;+
; NAME:
;       LNFIT
;
; PURPOSE:
;       LNFIT fits a log-normal distribution model to data.
;
; CATEGORY:
;       Data-Fitting
;
; CALLING SEQUENCE:
;       mu = LNFIT(x, sigma, XMIN=xmin, XMAX=xmax)
;
; INPUTS:
;       X -- The data to be fit.
;
; OPTIONAL INPUTS:
;       XMIN -- The pre-determined lower bound of the log-normal fit.
;       XMAX -- The pre-determined upper bound of the log-normal fit.
;
; KEYWORD PARAMETERS:
;       VERBOSE -- Be verbose in output.
;
; OUTPUTS:
;       MU    -- The mean of the log-normal distribution (i.e. ln[x])
;       SIGMA -- The width of the log-normal distribution (also in
;                ln[x])
;
; OPTIONAL OUTPUTS:
;       XMIN -- The algorithmically-determined lower bound of the
;               log-normal fit. 
;       XMAX -- The algorithmically-determined upper bound of the
;               log-normal fit. 
;
; COMMON BLOCKS:
;       LNFIT_BLOCK -- An internal COMMON block needed to pass
;                      information to the functions being
;                      minimized by Powell's method.
;
; NOTES:
;       This routine uses the POWELL function minimization routine,
;       which does not appear to contain anything in the way of
;       internal error handling.  Therefore, if POWELL reaches the
;       maximum number of iterations, it simply vomits all over the
;       session, essentially negating all computation up to that
;       point.  For stability, it seems keeping ftol at or below 1d-4
;       works -- not sure if step size is related to desired
;       tolerance.  If tolerance is too big, POWELL will step right
;       over the minimum and back without converging.
;
; MODIFICATION HISTORY:
;
;       Created:  04/23/13, TPEB -- Initial version, following
;                                   algorithm found in Olmi et
;                                   al. (2013, A&A, 551, A111).
;       Modified: 04/24/13, TPEB -- Routine is now functional, with
;                                   adequate documentation.
;       Modified: 05/08/13, TPEB -- Added functionality to specify a
;                                   hard upper and/or lower limit to
;                                   the range over which the
;                                   log-normal function is a valid
;                                   descriptor of the data.
;       Modified: 05/13/13, TPEB -- Fixed issue where POWELL must use
;                                   a COMMON block, and nested calls
;                                   to POWELL quickly run through the
;                                   available number of iterations.
;       Modified: 05/15/13, TPEB -- Have removed POWELL for the
;                                   fitting, since the IDL
;                                   implementation does not have any
;                                   error-chatching mechanism.  The
;                                   [mu,sig] fitting now handled by
;                                   TNMIN, with the [xinf,xsup] being
;                                   characterized by AMOEBA.  The
;                                   latter was constrained better to
;                                   not produce errors, and the former
;                                   utilizes built-in error catching
;                                   to avoid vomiting.
;       Modified: 09/04/13, TPEB -- Changed the sense in which AMOEBA
;                                   operates, in that the supplied
;                                   XMIN and XMAX are used as
;                                   *bounds*, but the routine may
;                                   select limits intermediate to the
;                                   provided values.
;
;-

;;===========================================================
;; Function to compute the exponential argument "x" of the log-normal
;;   distribution.  This argument crops up so much in the fitting of a
;;   log-normal distribution, that a separate function is most
;;   convenient.
FUNCTION LNFIT_X, M, mu, sig
  
  COMPILE_OPT IDL2, LOGICAL_PREDICATE, HIDDEN
  
  RETURN,(alog(double(M)) - double(mu)) / double(sig) / sqrt(2.d)
  
END


;;===========================================================
;; Function to compute the negative of the log-likelihood for a
;;   log-normal distribution.  This function is for use with TNMIN.
FUNCTION LNFIT_LOGL, arg, XFN=x, XINF=xinf, XSUP=xsup, $
                     _EXTRA=_extra
  
  COMPILE_OPT IDL2, LOGICAL_PREDICATE, HIDDEN
  
  ;; Parse out the argument into easily-usable variables.  
  mu   = double(arg[0])
  sig  = double(arg[1])
  xinf = double(xinf)
  xsup = double(xsup)
  
  ;; Normalization factor
  C = sqrt(2.d/!dpi) / sig / $
      ( erfc(LNFIT_X(xinf,mu,sig)) - erfc(LNFIT_X(xsup,mu,sig)) )
  
  ;; Log-Likelihood
  lnfitx = LNFIT_X(x,mu,sig)
  logL = ( double(n_elements(x)) * alog(C) ) - total( lnfitx * lnfitx , /DOUBLE)
  
  ;; For debugging purposes, print if values not finite
  IF ~finite(C) || ~finite(logL) THEN $
     print,string(mu,sig,xinf,xsup,C,logL,format="(F0.3)")
  
  ;; Return the log-likelihood for maximization
  RETURN,logL
  
END


;;===========================================================
;; Function for the exploration of [xinf,xsup] space with AMOEBA to
;;   find the set of values which minimizes the KS D-statistic
;;   between the data's CCDF and the model CCDF.  Given a value
;;   of [xinf,xsup], this function finds the maximum-likelihood
;;   [mu,sig] for a log-normal distribution using TNMIN, and computes
;;   the D-statistic between the data and that ML model.
FUNCTION LNFIT_LIM, arg
  
  COMPILE_OPT IDL2, LOGICAL_PREDICATE, HIDDEN
  
  COMMON LNFIT_BLOCK, xfn, fullx, xinf, xsup, ftol, fixcase
  
  ;; Parse out the argument into easily-usable variables, and in real
  ;;   space (rather than log space).
  CASE fixcase OF
     0: BEGIN                   ; No limits on (xmin,xmax)
        xmin = exp(double(arg[0]))
        xmax = exp(double(arg[1]))
     END
     1: BEGIN                   ; xmax bounded, but xmin is free
        xmin = exp(double(arg[0]))
        xmax = exp(double(arg[1])) < xsup
     END
     2: BEGIN                   ; xmin bounded, but xmax is free
        xmin = exp(double(arg[0])) > xinf
        xmax = exp(double(arg[1]))
     END
     3: BEGIN                   ; both (xmin,xmax) are bounded
        xmin = exp(double(arg[0])) > xinf
        xmax = exp(double(arg[1])) < xsup
     END
  ENDCASE
  
  
  ;; Find which data points fall within the passed limits.
  xind = where(fullx GE xmin AND fullx LE xmax, nx)
  IF ~nx THEN RETURN,!values.d_infinity
  xfn = double(fullx[xind])
  
  print,'XINF/XSUP',[xmin,xmax]
  
  ;; Set up to run TNMIN to determine the ML values of [mu,sig]
  P0 = [5.d,1.2d]
  parinfo = replicate( {value:0.d,$
                        fixed:0b,$
                        step:0.05d,$
                        limited:[1b,1b],$
                        limits:[double(alog((xmin/10.) < 0.1)),$
                                double(alog(xmax))]}, 2)
  parinfo.value = P0
  parinfo[1].limits  = [0.1d,100.d]
  fcnargs = {XFN:xfn, XINF:xmin, XSUP:xmax}
  
  P = TNMIN('LNFIT_LOGL', PARINFO=parinfo, /AUTODERIVATIVE, NITER=niter, $
            FUNCTARGS=fcnargs, NFEV=ncalls, STATUS=status, ERRMSG=errmsg,$
            QUIET=1, BESTMIN=bestmin, /ITERDERIV, /MAXIMIZE)
  IF status NE 1 && status NE -18 THEN print,status,' ',errmsg
  
  ;; Error-catching for TNMIN -- But -18 does not cause vomiting
  IF status LE 0 && status NE -18 THEN BEGIN
     message,'Catching TNMIN error; setting [mu,sig] = [5.0,1.2]',/inf     
     P = P0
  ENDIF
  
  ;; Replace values
  xinf = xmin
  xsup = xmax
  
  ;; Compute the data and model CCDFs.
  cx = ccdf(xfn)
  cf = ( erfc(LNFIT_X(xfn,P[0],P[1]))  - erfc(LNFIT_X(xsup,P[0],P[1])) ) / $
       ( erfc(LNFIT_X(xinf,P[0],P[1])) - erfc(LNFIT_X(xsup,P[0],P[1])) )
  
  print,[bestmin, max( abs(cf-cx) ),P,xinf,xsup]
  
  ;; Return the D-statistic -- which we want to minimize!
  RETURN, max( abs(cf-cx) )
END




;;**************************************************************************
;;==========================================================================
;; MAIN FUNCTION
;;**************************************************************************
FUNCTION LNFIT, x, sigma, XMIN=xmin, XMAX=xmax, VERBOSE=verbose, NAMOEBA=ncalls
  
  !except = 0                   ; Turn off annoying math warning messages!!!
  
  COMPILE_OPT IDL2, LOGICAL_PREDICATE
  
  COMMON LNFIT_BLOCK, xfn, fullx, xinf, xsup, ftol, fixcase
  
  verbose = KEYWORD_SET(verbose)
  ftol = 1.d-4
  
  ;; Keep a copy of the full distribution in the COMMON block
  fullx = x
  
  ;; Use AMOEBA to find limits
  fixmin  = n_elements(xmin) NE 0
  fixmax  = n_elements(xmax) NE 0
  fixcase = fixmax + 2*fixmin
  
  xinf = fixmin ? xmin : 1.0d1
  xsup = fixmax ? xmax : 1.d4
  
  ;; Choose case based on given limits (if any)
  message,'Using FIXCASE '+string(fixcase,format="(I0)"),/inf
  
  P0 = alog([xinf,xsup])
  nmaxiter = 100
  li = 0
  REPEAT BEGIN
     message,'Loop #'+string(++li,format="(I0)"),/inf
     print,P0
     P = AMOEBA( ftol, FUNCTION_NAME='LNFIT_LIM', FUNCTION_VALUE=fmin,$
                 NCALLS=ncalls, P0=P0, SCALE=2., NMAX=nmaxiter )
     nmaxiter += 100
     help,P,ncalls
  ENDREP UNTIL n_elements(P) EQ 2
  
  print,'P: ',n_elements(P)
  
  CASE fixcase OF
     0: BEGIN                   ; No limits on (xmin,xmax)
        xmin = exp(P[0])
        xmax = exp(P[1])
     END
     1: BEGIN                   ; xmax bounded, but xmin is free
        xmin = exp(P[0])
        xmax = xsup < exp(P[1])
     END
     2: BEGIN                   ; xmin bounded, but xmax is free
        xmin = xinf > exp(P[0])
        xmax = exp(P[1])
     END
     3: BEGIN                   ; both (xmin,xmax) are bounded
        xmin = xinf > exp(P[0])
        xmax = xsup < exp(P[1])
     END
  ENDCASE
  
  IF verbose THEN print,xmin,xmax,ncalls,fmin
  
  IF verbose THEN print,'LIMITS: ',[xmin,xmax]
  
  
  ;;===============================================================
  ;; Whether user-supplied or computed [xinf,xsup] values are used,
  ;;   now get ML value of [mu,sig].
  
  ;; Place limits into the COMMON block for use...
  xinf = xmin
  xsup = xmax
  
  ;; Find which data points fall within the passed limits.
  xind = where(fullx GE xinf AND fullx LE xsup, nx)
  IF ~nx THEN RETURN,!values.d_infinity
  xfn = double(fullx[xind])
  
  ;; Set up to run TNMIN to determine the ML values of [mu,sig]
  P0 = [5.d,1.2d]
  parinfo = replicate( {value:0.d,$
                        fixed:0b,$
                        step:0.05,$
                        limited:[1b,1b],$
                        limits:[double(alog((xinf/10.) < 0.1)),$
                                double(alog(xsup))]}, 2)
  parinfo.value = P0
  parinfo[1].limits[0] = 0.1d
  parinfo[1].limited   = [1b,0b]
  fcnargs = {XFN:xfn, XINF:xinf, XSUP:xsup}
  
  P = TNMIN('LNFIT_LOGL', PARINFO=parinfo, /AUTODERIVATIVE, NITER=niter, $
            FUNCTARGS=fcnargs, NFEV=ncalls, STATUS=status, ERRMSG=errmsg,$
            QUIET=0, BESTMIN=bestmin, /ITERDERIV, /MAXIMIZE)
  IF status NE 1 && status NE -18 THEN print,status,' ',errmsg
  
  IF verbose THEN print,'Log-Normal: ',P
  
  ;; Set function return values
  sigma = P[1]
  RETURN, P[0]
END
